//
//  Jacobian.cpp
//  EMT_C
//
//  Created by Duong Ngo on 8/27/16.
//  Copyright © 2016 Duong Ngo. All rights reserved.
//

#include "Jacobian.hpp"
#include "mapping.hpp"

#include <math.h>

void find_jacobian_structure_for_time_t(const int& t, Index* iRow, Index* jCol, int& dem)
{
    // Mapping f
    std::map<std::string, int> f;
    find_map_at_time_t(t, f);
    
    
    
    //Bankers
    //g0
    int pg= (t-1)*var_foc;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    
    //g1
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rf"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    
    //g2
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rm"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.mur"); ++dem;
    
    //g3
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.mur"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rn"); ++dem;

    //g4
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ql_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;

    //g5
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.n_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.tau"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.up"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.x"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.x_p"); ++dem;
    
    

    //g6
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rm_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.m_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.tau"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.n_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.up"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.x"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.x_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rn_p"); ++dem;

    //g7
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.mur"); ++dem;
    
    //g8
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    
    //g9
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;

   
    //Households
    
    //g10
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.etaz"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
 
    //g11
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    
    //g12
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rm"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lamb_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;

    //g13
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ql_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lamb_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.etab"); ++dem;

    //g14
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lamb_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pm_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.k"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.y_f"); ++dem;

    //g15
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rm_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.m_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.tau"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.etaz"); ++dem;
    
    
    //g16
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.etab"); ++dem;
    
    //g17
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pie_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.y_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pm"); ++dem;

    //g18
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.k_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.l"); ++dem;

    
    //g19
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.u"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pie_f"); ++dem;

  
    //g20
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;
    
    //g21
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.k"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.k_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    
    
    //g22
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;

    //g23
    ++pg;
    if (t< T_change_homotopy){
        iRow[dem]=pg; jCol[dem]=f.at("s.tau"); ++dem;
    }
    else {
        iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
        iRow[dem]=pg; jCol[dem]=f.at("s.m_p"); ++dem;
        iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;
    }
    
    //g24
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rn"); ++dem;

    
    //g25
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.l"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pm"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    
    //g26
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.up"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.up_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.w_f"); ++dem;
    
    //g27
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.x"); ++dem;
    
    //g28
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.w"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pm"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;

}

//***********************************************************************************
//***********************************************************************************

void find_jacobian_values_for_time_t (const int& t, const Number* x, Number* values, int& dem)
{
    // Var s
    variables s;
    find_variable_s_at_time_t(t, x, s);
    
    //Derivatives
    //Bankers
    //g0
    int pg= (t-1)*var_foc;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    values[dem]= s.cb;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    values[dem]= s.gamma;
    ++dem;
    
    //g1
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rf"); ++dem;
    values[dem]= -beta_b*s.gamma_f*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma_f"); ++dem;
    values[dem]= -beta_b*s.Rf*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    values[dem]= -beta_b*s.Rf*s.gamma_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    values[dem]= -2*pe*pow(fmax(s.muc,0),1);
    ++dem;
    
    //g2
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rm"); ++dem;
    values[dem]= -beta_b*s.gamma_f*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma_f"); ++dem;
    values[dem]= -beta_b*s.Rm*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    values[dem]= -beta_b*s.Rm*s.gamma_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    values[dem]= -2*pe*pow(fmax(s.muc,0),1);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.mur"); ++dem;
    values[dem]= -varphi*2*pow(fmax(s.mur,0),1);
    ++dem;
    
    //g3
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma_f"); ++dem;
    values[dem]= -beta_b*s.Rn*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    values[dem]= -beta_b*s.Rn*s.gamma_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    values[dem]= -2*pe*pow(fmax(s.muc,0),1);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.mur"); ++dem;
    values[dem]= -2*pow(fmax(s.mur,0),1);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rn"); ++dem;
    values[dem]= -beta_b*s.gamma_f*s.ip_f;
    ++dem;
    
    //g4
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    values[dem]= s.ql+thetaa;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    values[dem]= s.gamma;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ql_f"); ++dem;
    values[dem]= -beta_b*deltab*s.gamma_f*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma_f"); ++dem;
    values[dem]= -beta_b*(deltab+ deltab*s.ql_f)*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    values[dem]= -beta_b*(deltab+ deltab*s.ql_f)*s.gamma_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    values[dem]= -2*(1-vec_kappa[t])*pe*pow(fmax(s.muc,0),1);
    ++dem;
    
    //g5
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.n_p"); ++dem;
    values[dem]= s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    values[dem]= s.n_p;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.tau"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.up"); ++dem;
    values[dem]= s.x- s.x_p;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.x"); ++dem;
    values[dem]= s.up;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.x_p"); ++dem;
    values[dem]= -s.up;
    ++dem;
    
    //g6
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rm_p"); ++dem;
    values[dem]= -s.m_p*s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.m_p"); ++dem;
    values[dem]= -s.Rm_p*s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    values[dem]= -s.Rm_p*s.m_p+ deltab*s.bh_p+ (s.Rn_p-1)*s.n_p;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    values[dem]= -s.s;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;
    values[dem]= -s.ql;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    values[dem]= -thetaa;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh_p"); ++dem;
    values[dem]= deltab*s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.tau"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.n_p"); ++dem;
    values[dem]= (s.Rn_p-1)*s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.up"); ++dem;
    values[dem]= -(s.x-s.x_p);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.x"); ++dem;
    values[dem]= -s.up;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.x_p"); ++dem;
    values[dem]= s.up;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rn_p"); ++dem;
    values[dem]= s.n_p*s.ip;
    ++dem;

    
    //g7
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    values[dem]= -varphi;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.mur"); ++dem;
    values[dem]= 2*fmax(-s.mur,0);
    ++dem;

    //g8
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    values[dem]= (1-vec_kappa[t]);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    values[dem]= 1;
    ++dem;

    
    //g9
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh_p"); ++dem;
    values[dem]= -deltab*s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    values[dem]= -deltab*s.bh_p;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;
    values[dem]= -1;
    ++dem;
    
    
    //Households
    
    //g10
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.etaz"); ++dem;
    values[dem]= -2*pow(fmax(s.etaz,0),1);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    values[dem]= -1/s.ch/s.ch;
    ++dem;
    
    //g11
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    values[dem]= s.ch;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    values[dem]= s.lamb;
    ++dem;
    
    //g12
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rm"); ++dem;
    values[dem]= -beta_h*s.lamb_f*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lamb_f"); ++dem;
    values[dem]= -beta_h*s.Rm*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    values[dem]= -beta_h*s.Rm*s.lamb_f;
    ++dem;
    
    //g13
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    values[dem]= s.lamb;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    values[dem]= s.ql;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ql_f"); ++dem;
    values[dem]= -beta_h*deltab*s.lamb_f*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lamb_f"); ++dem;
    values[dem]= -beta_h*(deltab+ deltab*s.ql_f)*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    values[dem]= -beta_h*(deltab+ deltab*s.ql_f)*s.lamb_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.etab"); ++dem;
    values[dem]= -2*pe*pow(fmax(s.etab,0),1);
    ++dem;
    
    //g14
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lamb_f"); ++dem;
    values[dem]= -beta_h*(1-deltak);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pm_f"); ++dem;
    values[dem]= -beta_h*alphaa*s.lama_f*s.y_f/s.k;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama_f"); ++dem;
    values[dem]= -beta_h*alphaa*s.pm_f*s.y_f/s.k;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.k"); ++dem;
    values[dem]= -beta_h*alphaa*s.pm_f*s.lama_f*s.y_f*(-1)/s.k/s.k;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y_f"); ++dem;
    values[dem]= -beta_h*alphaa*s.pm_f*s.lama_f/s.k;
    ++dem;

    
    //g15
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rm_p"); ++dem;
    values[dem]= s.m_p*s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.m_p"); ++dem;
    values[dem]= s.Rm_p*s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    values[dem]= s.Rm_p*s.m_p- deltab*s.bh_p;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    values[dem]= s.s;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;
    values[dem]= s.ql;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.tau"); ++dem;
    values[dem]= 0;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh_p"); ++dem;
    values[dem]= -deltab*s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.etaz"); ++dem;
    values[dem]= 0;
    if (s.etaz < 0){
        values[dem]= -2*s.etaz;
    }
    ++dem;
    
    
    //g16
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.etab"); ++dem;
    values[dem]= 1;
    ++dem;

    
    //g17
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;
    values[dem]= -iota*2*s.pie+ iota;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama_f"); ++dem;
    values[dem]= iota*beta_h/s.lama*(s.pie_f-1)*s.pie_f*s.y_f/s.y;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    values[dem]= iota*beta_h*s.lama_f*(-1)/s.lama/s.lama*(s.pie_f-1)*s.pie_f*s.y_f/s.y;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pie_f"); ++dem;
    values[dem]= iota*beta_h*s.lama_f/s.lama*(2*s.pie_f-1)*s.y_f/s.y;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y_f"); ++dem;
    values[dem]= iota*beta_h*s.lama_f/s.lama*(s.pie_f-1)*s.pie_f/s.y;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    values[dem]= iota*beta_h*s.lama_f/s.lama*(s.pie_f-1)*s.pie_f*s.y_f*(-1)/s.y/s.y;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pm"); ++dem;
    values[dem]= epsilon;
    ++dem;
    
    //g18
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.k_p"); ++dem;
    values[dem]= -alphaa*pow(s.k_p, alphaa-1)*pow(s.l,1-alphaa);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.l"); ++dem;
    values[dem]= -(1-alphaa)*pow(s.l,-alphaa)*pow(s.k_p,alphaa);
    ++dem;
    
    
    //g19
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.u"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pie_f"); ++dem;
    values[dem]= -1/beta_b*phipie*pow(s.pie_f, phipie-1);
    ++dem;
    
    
    //g20
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    values[dem]= 1- iota/2*(s.pie-1)*(s.pie-1);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    values[dem]= -thetaa;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;
    values[dem]= -iota*s.y*(s.pie-1);
    ++dem;
    
    //g21
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.k"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.k_p"); ++dem;
    values[dem]= -(1-deltak);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    values[dem]= -1;
    ++dem;
    
    
    //g22
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    values[dem]= s.pie;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;
    values[dem]= s.ip;
    ++dem;
    
    //g23
    ++pg;
    if (t< T_change_homotopy){
//        iRow[dem]=pg; jCol[dem]=f.at("s.tau"); ++dem;
        values[dem]= 1;
        ++dem;
    }
    else {
//        iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
        values[dem]= 1/s.m;
        ++dem;
//        iRow[dem]=pg; jCol[dem]=f.at("s.m_p"); ++dem;
        values[dem]= -1/s.m_p;
        ++dem;
//        iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;
        values[dem]= (1+rho_m)/s.pie;
        ++dem;
    }


    
    //g24
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rn"); ++dem;
    values[dem]= 1;
    ++dem;
    
    //g25
//    iRow[dem]=pg; jCol[dem]=f.at("s.l"); ++dem;
    values[dem]= chi*(nu+1)*pow(s.l,nu);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pm"); ++dem;
    values[dem]= -(1-alphaa)*s.y*s.lama;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    values[dem]= -(1-alphaa)*s.pm*s.lama;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    values[dem]= -(1-alphaa)*s.pm*s.y;
    ++dem;
    
    //g26
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    values[dem]= s.up;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.up"); ++dem;
    values[dem]= s.lama;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama_f"); ++dem;
    values[dem]= -beta_h*(s.up_f+ s.w_f);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.up_f"); ++dem;
    values[dem]= -beta_h*s.lama_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.w_f"); ++dem;
    values[dem]= -beta_h*s.lama_f;
    ++dem;
    
    //g27
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.x"); ++dem;
    values[dem]= 1;
    ++dem;
    
    //g28
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.w"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pm"); ++dem;
    values[dem]= s.y;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;
    values[dem]= iota*s.y*(s.pie-1);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    values[dem]= -(1-s.pm)+ iota/2*(s.pie-1)*(s.pie-1);
    ++dem;

    
    
}

//***********************************************************************************
//***********************************************************************************


void find_jac_struct_for_initial_condition (Index* iRow, Index* jCol, int& dem)
{
    // Starting position
    int pg = (T-1)*var_foc; //Staring of equation
    int px = 0; // Staring variable;
    
    for (int i=0; i<var_foc; ++i){
        iRow[dem]=pg+i; jCol[dem]=px+i; ++dem;
    }
}

void find_jac_values_for_initial_condition (const Number* x, Number* values, int& dem)
{
    for (int i=0; i<var_foc; ++i){
//        iRow[dem]=pg+i; jCol[dem]=px+i; ++dem;
        values[dem]= 1;
        ++dem;
    }
}

//***********************************************************************************
//***********************************************************************************

void find_jac_struct_for_terminal_condition (Index* iRow, Index* jCol, int& dem)
{
    // Starting position
    int pg = T*var_foc; //Staring of equation
    int px = T*var_foc; // Staring variable;
    
    for (int i=0; i<var_foc; ++i){
        iRow[dem]=pg+i; jCol[dem]=px+i; ++dem;
    }
}

void find_jac_values_for_terminal_condition (const Number* x, Number* values, int& dem)
{
    for (int i=0; i<var_foc; ++i){
//        iRow[dem]=pg+i; jCol[dem]=px+i; ++dem;
        values[dem] = 1;
        ++dem;
    }

}

//***********************************************************************************
//***********************************************************************************

void summarize_jac_struct (Index* iRow, Index* jCol, int& dem)
{
    for (int t=1; t<T; ++t){
        find_jacobian_structure_for_time_t(t, iRow, jCol, dem);
    }
    find_jac_struct_for_initial_condition(iRow, jCol, dem);
    find_jac_struct_for_terminal_condition(iRow, jCol, dem);
}

void summarize_jac_values (const Number* x, Number* values, int& dem)
{
    for (int t=1; t<T; ++t){
        find_jacobian_values_for_time_t(t, x, values, dem);
    }
    find_jac_values_for_initial_condition(x, values, dem);
    find_jac_values_for_terminal_condition(x, values, dem);
}






